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A general numerical method is presented to locate the partition function zeros in the 
complex j3 plane for large lattice sizes. We apply this method to the 2D Ising model and 
results are reported for square lattice sizes up to L = 64. We also propose an alternative 
method to evaluate corrections to scaling which relies only on the leading zeros. This 
method is illustrated with our data. 



1. Introduction 

Although the two-dimensional Ising model is exactly solved for zero external field, 
it continues receiving attention in many aspects. In a recent work by Beale 1 for 
instance, the low-temperature series expansion for the partition function was ex- 
actly determined for finite lattices with periodic boundary conditions. In terms 
of the expansion variable u — e~ 4/3 the partition function on a mxn lattice size 
becomes a polynomial of finite degree in u, and its coefficients g(E), the number of 
configurations with energy E, were calculated from an exact closed form based on 
Kaufman's solution. 2 

Motivated by that calculation as well by the enhancement of computer facilities 
we decided to revisit the 2D Ising model to obtain the exact partition function 
zeros in the complex temperature plane. 3 This approach had already been pursued 
by Katsura and Abe 4 ' 5 in the early investigations of zero distributions in order to 
check the proposal by Fisher about their loci. 6 Other papers have aimed the study 
of the critical properties from the leading zeros, but all of them were limited to 
small lattices (mxn < 13x13). 7,8,9 

In this work we present a procedure to obtain exact complex zeros of the partition 
function for large lattice sizes. We provide a description of a way round technical 
limitations on solving polynomials, at least in what concerns the location of the 
first zeros (it?, wP>, ■■•)■ It is a modified version of the scanning procedure which 
has been applied to continuous energy distributions of lattice gauge theories. 10,11,12 
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Our approach is based on constructing a new function in terms of \ng{E). It is 
presented in section 2, where the first exact zeros are calculated for square lattice 
sizes up to L = 64. Since we are working with numerical computation, we mean by 
exact zeros accurate values limited only by the use of double-precision floating-point 
arithmetic. 

The partition function zeros approach has been largely used to obtain infor- 
mation on phase transitions from Monte Carlo (MC) simulations 3,10,13 or exact 
enumerations 9 of finite systems. In this context, Itzykson's et al. 3 finite size scaling 
(FSS) relation for the first zero is the proper way to calculate the correlation expo- 
nent v. We present results so obtained in section 3. In addition, looking forward to 
obtain more information from those exact zeros, we propose, in the same section, a 
new way to evaluate corrections to FSS which relies only on u\(L) data. 

2. Exact partition function zeros 

The partition function of the two-dimensional Ising model on a m x n lattice can be 
written as a polynomial, 



where u = e 4/3 . 

Kaufman's solution for the isotropic Ising model renders the analytical expres- 
sions to be expressed in the polynomial form (|l|). Following Beale 1 , this can be 
done for any lattice size by using MATHEMATICA. 

A further step namely exact determination of their zeros and FSS analysis for 
the leading ones, can be achieved with this polynomial form. However we have 
checked that it was not possible to handle systems for L larger than 16 with our 
workstations. In fact, as the lattice size increases, the exact coefficients become 
very large integers. The enormous increasing of their maximum values, typically 
\ng(E) ~ 174 for L = 16 and going up to \ng(E) ~ 2835 for L — 64, prevent us 
from solving them by using computer algebra language. The same reasons do not 
indicate the use of standard numerical algorithms, 14,15 usually employed in those 
cases. 

To circumvent this problem we borrowed inspiration from lattice gauge theories 
where the energy distributions are continuous. There, in contrast to spin systems 
where the action takes discrete values and the partition function becomes a poly- 
nomial in u, a time series analysis in function of the complex coupling /? = f3 x + i(3 y 
is more efficient in calculating the first zero which is closest to the infinite volume 
critical point. It is a two step approach. 12 First, we scan graphically for separate 
zeros of ReZ((3) and lmZ((3), where 



nin 
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The shift AE = E — < E > is usually introduced for technical reasons, to 
avoid numerical overflow computations, although it is not any more relevant in 
the new approach. A typical output is shown in Fig. 1. Crosses correspond to 
Re Z{(5) = and diamonds to Im Z{j3) — 0. The wanted zero is obtained when 
the lines cross. Second, we compute this zero to a desired precision as an iterative 
process. This can be achieved by means of the minimization algorithm AMOEBA 16 
for the function ((Re Z{f3)) 2 4 (Im Z(/3)) 2 ) 1 / 2 , whose starting point is obtained 
from a simple inspection of figures like Fig. 1. As an example, we can use the input 
(0.43765, 0.0131) as the starting point to this routine which leads, after roughly 100 
iterations, to (/?°, /3°) = (0.4376431265, 0.01311604331). 

Now we shall describe how to implement our approach. 

Since our aim is to achieve large lattices one has to work with logarithms. For 
this end we need to introduce a new function F{(3) — F x ([3) 4 iF y {(3) to play the 
role of Z((3) itself. We start from splitting Re Z(f3) into two positive parts, namely 
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Fig. 1. Search for the first partition function zero for L = 64. The crosses indicate the zeros for 
ReF(/3) and the diamonds the ones for ImF(/3). The complex function F(0) has the same zeros 
as Z(j3). 



Ge and He, defined by 



Re Z([3) — Ge — He 



(3) 
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where 

G E = ' 9(E) e-^ AE cos(a E ) , (4) 

E 

H E = Y / " 9(E) e-^ AE \cos{a E )\ , (5) 

E 

and a E = A(3 y /\E. Here means a summation over E provided cos(a£) > 0, 
and J2" stands for the complementary values where cos(o:_e) < 0. Next wc 
calculate lnG_E and lnH E in a recursive way from the terms In g(E) , lncos(a_E) 
and — 4(} X AE. Logarithmic terms can be added up two by two in and H E , 
respectively, by using the relation 

In(o + 6) = ln& + ln(l + e lna - lnb ). (6) 

Finally, since we are interested in the (3 values where ReZ(f3) and ImZ(p) change 
signals, we realize that this can be achieved by the function 

F x {(3)= In G E — In H E , (7) 

and a similar one for F y ((3), which follows from the imaginary part of Z(f3). 

Now, we apply the two steps procedure to find roots of F x {(3) and F y {f}) instead 
of ReZ{f3) and lmZ{/3). In Fig. 1 wc show the first step for L = 64. Table 1 
contains our leading zeros u\{L) for lattice sizes up to L = 64, with rounded errors 
in the last digit. 



Table 1. First partition function zeros. 



L 


Rc(«») 


Im(u^) 


4 


0.1624473772 


0.16648190032 


6 


0.1756913616 


0.10528348725 


8 


0.1780809275 


0.07710375572 


9 


0.1783370200 


0.06801661701 


10 


0.1783571854 


0.06084948478 


12 


0.1780873239 


0.05026266796 


15 


0.1774653671 


0.03986421638 


16 


0.1772557409 


0.03729300267 


18 


0.1768587209 


0.03303236187 


20 


0.1764984476 


0.02964570468 


24 


0.1758873488 


0.02460155919 


30 


0.1751918649 


0.01959967784 


32 


0.1750048586 


0.01835571819 


36 


0.1746815338 


0.01628818111 


40 


0.1744124041 


0.01463927731 


48 


0.1739912845 


0.01217440355 


60 


0.1735492820 


0.00971963267 


64 


0.1734355215 


0.00910750889 


oo 


0.1715728753 





This method can easily be implemented and takes few minutes of CPU time 
for a Fortran code in a workstation after we have calculated the coefficients g(E). 1 




The lattice sizes were chosen to explore finite size corrections by an alternative 
method which will be presented in the next section. In addition, we plot in Fig. 2 
all zeros for L = 16 obtained with MATHEMATICA and compare them with the 
expected phase boundaries in the u plane for L — > oo. 17 ' 18 This curve corresponds 
to the locus of points where the free energy is non-analytic and it is parametrized by 
Re(u) = l+2 3 / 2 cosw + 2cos2w and Im(u) = 2 3 / 2 sin w + 2 sin2w, for < w < 2ir. 

3. Finite size scaling analysis 

The systematic dependence of v on finite systems can be explored to evaluate the 
main correction to scaling. From pairs of lattices L and L', we define the corre- 
sponding finite size estimators, 

This equation was already used to estimate the critical exponent v for the 3D Ising 
model from increasing pairs of lattices. 19 

In Table 2 we present sequences of 1/vl,sL as a function of a fixed rescaling 
factor s = L'/L. As min(L, L') increases, the values obtained by matching pairs of 
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Table 2. Sequence of estimates for 1/vl,sL- 



L / s 


1.5 


2 


3 


4 


4 


1.131948204 


1.107535391 


1.083907993 


1.071996844 


6 


1.067288812 


1.055806270 


1.044309560 


1.038354523 


8 


1.043516688 


1.036458296 


1.029248898 


1.025453193 


10 


1.031641084 


1.026690577 


1.021565572 


1.018835580 


12 


1.024655805 


1.020902776 


1.016977082 


1.014867660 


16 


1.016924350 


1.014448089 


1.011818989 


1.010388613 


20 


1.012804317 


1.010980582 


1.009024731 




24 


1.010266064 


1.008832544 






32 


1.007324514 


1.006329138 






40 


1.005681181 









lattices approach the expected limiting value v = 1. 

Equation (||), beyond being quite important to estimate v can be a starting point 
to evaluate finite size corrections, which can be due to a variety of sources. 20,21 For 
this end let us briefly recall Nightingale's finite size RG transformation. 

Under the hypothesis the system is large enough to consider the scaling rela- 
tion for the longitudinal correlation length £l(/?)> the standard expression for the 
correlation exponent v is 22,23 



where the scaling equation for the correlation length is given by 

a = LY t ((/3 - p^L 1 '" \hW* \uV*) . (10) 

This differentiable equation includes corrections due to the leading bulk irrelevant 
scaling field u with exponent j/3 < 0, and a magnetic field dependence for the sake 
of completeness. 

From Eq. (§) and Eq. © one obtains, for h = 0, 

f f L' V3 - V> 3 L' 2y3 - l 2 ^ , s 

a : 77777T + b — ——— + ... (11) 



vl,v v \n{L'/L) \n(L'/L) 

where ao and 60 include derivatives like dY% (y, z)/dy\ y =o tZ =o- If we replace L 
L' — 1 in the above equation, one obtains Privman and Fisher's results. 23 

Our analysis follows with the introduction of the rescaling factor s in Eq. (|Tl 
The exponent v can be asymptotically obtained from sequences either by extrap- 
olating In s — * 00 or by extrapolating In s — > (see the similar scaling behaviour 
of Binder's function 24 W% for obtaining 2j3/v). The exponent 2/3 can be evaluated 
through a linear regression for finite lattice sizes L if we fix the ratio s. At this 
point we call attention to the known fact that the main finite size correction for the 
2D Ising model comes from nonlinear scaling fields 23 ' 25 and gives origin to different 
L— dependent corrections in Eq. (pr|). 
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Table 3. Estimates of w from 1/i>l,sL data in Table 2. w refers to 1/vl.sL data obtained with 
the replacement \u®(sL) — u c \ / \u1(L) — by Im u®(sL) / Irau^L) in Eq. 8. 



s L —w — lna(s) — w — lna(s) 



1.5 


(4,6,8,10,12,16,20,24,32,40) 


1. 


.34903 


0.284605 


1 


.09706 


0.586739 


1.5 


(12,16,20) 


1. 


.28388 


0.514528 


1 


.05752 


0.716982 


1.5 


(12,16,20,24,32,40) 


1. 


.21665 


0.698876 


1 


.03807 


0.770224 


1.5 


(20,24,32,40) 


1. 


.17201 


0.850741 


1 


.02581 


0.811914 


1.5 


(32,40) 


1. 


.13858 


0.970497 


1 


.01865 


0.837543 


2 


(4,6,8,10,12,16,20,24,32) 


1. 


.35083 


0.458201 


1 


.09585 


0.735457 


2 


(12,16,20) 


1. 


.26133 


0.735541 


1 


.05021 


0.873225 


2 


(12,16,20,24,32) 


1. 


.21780 


0.853411 


1 


.03808 


0.906048 


2 


(20,24,32) 


1. 


.17102 


1.00516 


1 


.02555 


0.946659 


3 


(4,6,8,10,12,16,20) 


1. 


.37931 


0.624774 


1 


.10630 


0.903252 


3 


(12,16,20) 


1. 


.23804 


1.00130 


1 


.04348 


1.06970 


1 


(4,6,8,10,12,16) 


1. 


.39591 


0.736397 


1 


.11314 


1.01358 


4 


(12,16) 


1 


.24609 


1.11215 


1 


.04606 


1.18084 



Let us call w the effective exponent coming from the equation 

In f — ) =wlnL + lna(s), (12) 

\VL,sL V) 

which intends to detect the main correction regardless its origin. 

We collect in Table 3 our results for w. The third and fourth columns correspond 
to fit Eq. © to data of Table 2. 

For the 2D Ising model u c is exactly known, however for many models the value 
of u c is not known with high precision and in this case it is usual to replace \u\ — u c \ 
by Im u\. For sake of illustrative purposes both cases were considered. In our tables 
we use the notation w instead of w when the 1/vl.sL data is obtained from Eq. (^) 
with the replacement \u\(sL) — u c \ / \u°(L) — u c \ by Im u1(sL) / hxiu^(L), Different 
fixed ratios s are used to show the behaviour of w. As s increases, and consequently 
i', the numerical results show a trend in direction of w = — 1. Large s means 
working with crossings involving a large L, hence close to C . The corresponding 
best linear fits for all data are presented in Fig. 3. It is clear that small lattice sizes 
give origin to deviations in the employed linear equation (|l2|). 

To complete our analysis we present in Table 4 the dependence of w on small 
lattices. The value w w —1.7 is quite close to Binder's reported value, w ms —1.8, 
as the main correction to the function W%. 24,26 However there smaller lattices were 
considered in the analysis which seem to increase \w\. This trend can be caught 
from our smaller data set with s = 1.5 for L = 4 and 6. 

In summary, we have described an approach to compute partition function zeros 
for large lattice sizes. Although it was applied to the exact 2D Ising model partition 
function, it can be quite useful either when we deal with large coefficients, even out of 
scope of double precision computations, or when the polynomial has a large number 
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Fig. 3. Linear regressions for In {1/i^l,sL ~ 1) f° r several values of the rescaling factor s, according 
to Eq. (12). 

of coefficients, which would prevent us from using standard solving algorithms. This 
last situation is proper for the multicanonical simulation of the density of states in 
the 3D Ising model. 27 

In MC simulations FSS behaviour of the first zero has been used to estimate 
the exponent v and the critical coupling fi c . Beyond any MC data, with limited 
statistical precision, we were able to explore here the performance of the FSS anal- 
ysis proposed to evaluate corrections to scaling. This approach reveals to be quite 



Table 4. Estimates of w as in Table 3. 



s 


L 


— w 


—In a(s) 


—w 


— lna(s) 


1.5 


(4,6) 


1.66085 


-0.277077 


1.27701 


0.268877 


1.5 


(4,6,8) 


1.60430 


-0.191966 


1.23731 


0.328639 


1.5 


(4,6,8,10) 


1.56144 


-0.123144 


1.20877 


0.374454 


2 


(4,6) 


1.61773 


-1.27172 xl0~ 2 


1.24410 


0.478115 


2 


(4,6,8) 


1.56425 


6.77855 xl0~ 2 


1.20847 


0.531757 


2 


(4,6,8,10) 


1.52369 


0.132909 


1.18311 


0.572470 



useful in setting an upper limit on w for the same data set as s increases, as can be 
observed, for example, from our Table 3 for L = 12, 16 and 20. 
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